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TECHNICAL MEMORANDUM 


VARIATIONAL DIFFERENTIAL EQUATIONS FOR ENGINEERING 
TYPE TRAJECTORIES CLOSE TO A PLANET WITH AN ATMOSPHERE 


INTRODUCTION 

This report documents the theory behind a subroutine for optimal 
trajectory computation. The mechanics involved are considered to be standard 
knowledge in atmospheric flight mechanics; therefore, no derivation of the 
equations of motion will be given. The subroutine is intended for use in 
connection with extremal field methods. Since a large amount of computational 
work load stems from the nonalignment of the inertial and the aerodynamic 
velocity vector, a special path in the subroutine has been provided for crude 
first investigations with aligned velocity vectors (no planet rotation, no wind) . 

The subroutine is intended for engineering -type applications. This 
determines the degree of sophistication for the mathematical model which is 
described in the next section (Mathematical Model) , In the computation of 
optimal trajectories, the selection of coordinate systems in which the motion 
is to be described is even more crucial than it is for straightforward numerical 
integration. It affects not only the computational work load (e. g. , use of 
trigonometric functions) and the integration accuracy but also the range of 
validity of linear approximations used in the iteration process and the interpret- 
ability of the adjoint variables. The coordinate systems chosen are described 
in the third section (coordinate Systems) . 

Symbols chosen in the report agree, in general, with those coded. The 
code is intended to be mnemonic (which, of course, is strongly biased by the 
author's background). 


MATHEMATICAL MODEL 


Experience with reentry trajectory computations [ 1, 2, 3] led to the 
definition of the following mathematical model. 



Planetary Model 


Only the gravitational force of the central body is considered. The J2- 
gravitational potential spherical harmonic term is taken into account in the 
vertical gravity component. Other influences of the planets’ gravitational field 
or of moons, other planets, and the sun are disregarded. 

The geometrical form of the planet is a flattened sphere. The actual 
shape is approximated by a third-order polynomial fit to the curve Ar 0 (A) * 
radial deviation from a sphere as function of the geocentric latitude for 
0 s A s i/2, Its derivative for A = 0 and n/2 is 0. For the earth the 
radius difference between pole and equator is approximately 21 km with a 
maximum in the slope between sphere and ellipsoid of almost 0. 2 deg around 
A = 7r/4. Therefore, for shallow entry angles this has to be taken into account 
12 ]. 


Atmospheric characteristics are a function of altitude alone. Density, 
pressure, and velocity -of-sound profiles have to be furnished by subroutines 
which also provide their derivatives. These profiles are considered to be 
independent of longitude, latitude, and time in the ascent or reentry region 
traversed. 


Vehicle Model 

The vehicle is considered to be a point mass. It has no moments of 
inertia, and its angular motions are considered to be instantaneously change- 
able controls. 

There acts no aerodynamic side force. The vehicle is considered to 
have a plane of symmetry in which the aerodynamic force vector is bound to 
lie (no sideslip). 

The aerodynamic characteristics of the vehicle are to be furnished by 

a subroutine. The program is set up for a drag polar type of representation 

(C = C + CjC + C 2 C£. + C 3 C?_). Provision is made to take into account 
JD Oq Ju jL JL 

Mach-number-and altitude- (viscous interaction) dependent changes of the 
characteristics. Note that this does not allow for double-valued polars (back- 
side of the lift-drag curve) . This case may be dealt with by taking the drag 
coefficient as control instead of the lift coefficient or by going to a parameter 
representation (angle of attack) , which, however, doubles the amount of 
computation necessary. 
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Thrust is provided by rocket engines with constant mass flow. The 
effect of atmospheric back pressure is taken into account. 

The thrust angle as control is considered to be unconstrained and not 
linked to a certain body axis. This yields information which might be taken into 
account for the design of the vehicle. Thrust angle constraints may be easily 
added. 


COORDINATE SYSTEMS 

The basic plane of reference is the equatorial plane of the planet. 

Inertial Coordinate System (Index i) 

A Cartesian planet-centered coordinate system with fixed inertial 
directions is considered to be an inertial system, neglecting the influence of 
the planets motion around the sun or largerscale motions. The x - and y.- 

axes lie in the equatorial plane, and the planet rotates around the z. axis 

(Fig. l) . The vernal equinox is taken as the reference direction for x^. 



Figure 1. Inertial and planet fixed coordinate systems. 
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Planet-Fixed Coordinate System 

This system, indexed p, rotates relative to the inertial system around 
the common z-axis. Its reference direction is taken to be the zero-meridian 
crossing with the equator. In this system the position of the vehicle is measured 
in polar coordinates: longitude 0, latitude A, and radius r. The altitude 
above sea level is computed from r and A . The rate of change of altitude is 
not equal to r as it is for a spherical planet. 


Radius Normal Coordinate System (Index g) 

The equations of motion are written in a Cartesian coordinate system 
the z-axis of which is aligned with the radius vector from the planet mass- 
center to the vehicle center of gravity. The x-axis is in the direction of the 
inertial velocity component normal to the radius vector (Fig. 2) . By this 
definition, there is no velocity component in the y-direction. Forces normal 
to the momentary plane of motion containing the planet mass center turn the 
’’horizontal" component U of the velocity vector and thereby the orientation of 
the x-axis. The direction of the x-axis is measured by the inertial azimuth 
angle x. from north positive over east. 

This choice of coordinate system is a compromise with respect to the 
factors mentioned in the introduction. For the most common reentry and the 
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final part of ascent trajectories, the kinetic energy of the vehicle is almost 
entirely represented in the horizontal velocity component U (small flight path 
angles) . The term U is an essentially monotonic function compared to oscil- 
latory functions for the velocity components in a Cartesian inertial or longitude- 
latitude oriented formulation. However, the azimuth angle is oscillatory and 
invokes trigonometric functions. There is a singularity in crossing the poles 
resulting in A changing sign and 0 jumping by it. The important quantity 
radial velocity is a state variable directly. 


Aerodynamic Coordinate System (Index a) 

The aerodynamic force coefficients are usually defined in a coordinate 
system, the x-axis of which is aligned with the aerodynamic velocity vector. 
The z-axis lies in the vehicle plane of symmetry. The drag cooefficient 

acts in the negative x - direction, the lift coefficient C is here defined 

3 ij 

upward for p = 0, where j u is the bank angle between the z-axis and the 
vertical plane. The direction of the x^ - axis with respect to the radius 

normal coordinate system is given by the aerodynamic flight path angle y in 

3 

the vertical plane and the horizontal misalignment angle in the azimuth Ay 

3 

between the horizontal inertial and aerodynamic velocity components (Fig. 3) . 


EQUATIONS OF MOTION AND HAMILTONIAN FUNCTION 

In this section the differential equations for the state variables will be 
given together with transformation relations. 



Figure 3. Aerodynamic coordinate system. 


5 






Transformation Aerodynamic -Radius Normal 

From Figure 3 it is seen that the north-south and east-west aerodynamic 
velocity components are given by 


V = U cos X + W 
NS 1 NS 


V EW = U Sln X i + W EW - r ^p COS A • 


( 1 ) 

( 2 ) 


where the wind components W are positive for north wind and east wind. The 
total horizontal aerodynamic velocity component is 


V ah= J V ki + V W“| U!+2U[ c° s *i 


w. 


NS 


+ Sin X i (W EW “ r “p C ° S A)] + W NS 
+ (W EW “ r % COS A)2 | ^ ‘ 


(3) 


The vertical aerodynamic velocity component with W positive for a downwind 
is 


V = r + W 
v v 


(4) 


yielding a total velocity vector 


V = (v 2 + v 2 + V 2 ) 

a v NS EW v ; 


% 


(5) 


The case of vertical winds will not be pursued further (W = 0) . The 
areodynamic flight path angle is given by 


sin y = V /V 
a v a 


or tan y = V /V , 
a v ah 


( 6 ) 


From Figure 3 one reads 



which by trigonometric relations yields 


sin Ax a = (V EW cos X, - V Ng sin X,)/V ah = F(U. X,, r, A) 

cos Ax a = (V EW sin X. + V Ng cos X^/V^ = F(U, Xj. r, A) (7) 


for the horizontal misalignment angle Ay . These expressions depend on the 

3 

state variables U, X.» r, A as indicated after the second eqvial sign. 

Forces are to be transformed from the aerodynamic into the radius 
normal coordinate system by 


X 

g 


’ll 




Y 

g 

= 

h 



X " 

a 

Z 

g 


*3 



z 

a 


where 


( 8 ) 


L = cos y cos Ay , 
a a 

U = cos y sin Ax , 
a a 

1 3 = -sin y , 

3 

m = cos m sin y cos Ax + sin M sin Ax , 

1 a a a 

= cos i u sin y sin Ax - sin n cos Ax , 

3 3 3 

n 3 = cos ix cos y 

3 

The bank angle ju is a control to be determined such that a pay-off function is 
extremized. 
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State Differential Equations 

Aerodynamic and thrust terms 1 are grouped together: 

* * Ur T 

Xi = U = - — + — cos e cos a 


r m 


- — [C cos y cos Ax 
_m_ D a_ a 

+ C (cos M sin y cos Ax + sin n sin Ax )] , ( 9 ) 

L H 3 3 

* * IT A T 

X, = x .= — sin x. tan A + — - sin e 
* i r l mU 

[CX cos y sin Ax 
Um __D a a 

+ C T (Cos n sin y sin Ax - sin / x cos Ax )] , (10) 

Xj 3 3 3 


li 

CO 

• « 

r " 

U 2 

r 

GM 

1 _ j(^-) 2 (3 sin 2 A - l) 

^ T 

+ — cos e sin c r 
m 



_Sq 

m 

‘ C _D 

sin y - C_ cos y cos p] 
a L a 

, (11) 

X 4 = 

r = 

*2 

9 


(12) 

X 5 = 

A = 

u 

— cos X. 
r i 

9 

(13) 

X S = 

m = 

= -b 

9 


(14) 



U sin x. 



x 7 = 

$ = 

1 

r cos A 

- C 0 

p 

(15) 


Aerodynamic terms will be designated by broken underlining, thrust 
terms by double underlining, and inertial terms by single underlining. 
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The lift coefficient is the second control. The drag coefficient is 

dependent on Mach-Reynolds number and the lift coefficient. A representation 
of the form 

C D = °D„ + C > C L + °’ C ’l + C »° 3 L (16) 

is implemented, where C and the C. contain the influence of the Mach number 

D 0 1 

or the combined Mach-Reynolds number effects (viscous interaction) . This 
form is adapted to a bivariate spline-representation of the drag polar. The 
term S is the aerodynamic reference area. The dynamic pressure 

q - \ p v* a (17) 

contains five state variables [see eqs. (l) through (5)] which are U, r, r, 

X.j A , where r and A also determine the air density p (h) . The altitude h 

above the planet reference surface is written 

h = r - R 0 + a. A , 1 = 1, 3 (18) 

with R 0 = equivalent spherical planet radius (6371.2 km for earth). The 
thrust T contains an atmospheric term 

T = b ' u e - (p « ' V ’ <19) 

where 

U e = exit velocity , 

P ro = free-stream pressure , 

A^ = exit area of engines 

It is decomposed into a component in the vertical plane containing X (thrust 

S 

angle cr positive from horizontal upward) and a component normal to this 
plane (angle e, positive for increasing x^) . 
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The Hamiltonian Function 


The adjoint variables are symbolized by the state variable symbol as 
index (rp for r) . The abbreviations used in the derivation are also used as 
code in the FORTRAN program. The Hamiltonian is decomposed into inertial, 
thrust, and atmospheric terms: 


H ■ H I + h t + h a 


Inertial terms: 


Hi - ~ (H4) - * rp ^ [l - *(2?) 2 (3 


sin' 


L 2 A - l) 


] 


+ X r - X.o? 
r 9 ,p 


( 20 ) 


( 21 ) 


where 


H 4 =-(HI2) + sin x. 


+ X. cos x. 
A i 


HI2 = X r - X U , 
u rp 


HI5 = X tan A + X„/cos A 
X 0 


Thrust terms: 


/ U 

H_ = b ( — w - X 
T V m m 


) 


A 

*oo e 
m 


w 


( 22 ) 


where 


w - Wi cos e + x . 

1 ~ sin e 
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Wi = X cos cr + X sin cr 
1 u rp 

In order to determine the optimal thrust program, has to be extremized 

with respect to the mass flow b and the thrust angles cr and e. The last 
term in equation (22) is dependent on the atmosphere too. 
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Atmospheric terms: 

P OT A „ 

H a = w — — ^ [C HI + C (H2 cos ju + H 3 sin m)1 , (23) 

A m m D L 

where 

HI = HO cos 7 + A sin y , 

a rp a 

H2 = HO sin 7 - A cos 7 , 

a rp a 

A 

H3 = A sin Ax — cos Ax , 

u si U a 

and 

HO = A cos Ax + A /U • sin Ax 
u a x a 

The term in square brackets is symbolized by 

H 6 = C HI + C T (H2 cos M + H3 sin n) . (24) 

13 Li 

The extremization of H 6 with respect to the controls C T (or C_, 01 ) and 

Jj D 

the bank angle n yields candidates for the optimal controls. The first term in 
equation (23) is an atmosphere-dependent thrust term. 

Reduced Equations for Aligned Inertial 
and Aerodynamic Velocity Vectors 

The misalignment angle Ax is identically zero, and aerodynamic 

a 

and inertial flight path angles are identical: 

7 = 7 .= tan - 1 (r/u) . (25) 

a 1 

With 

V - (r 2 + U 2 ) 1/2 , (26) 


11 



there is 


sin y = r/V. ; cos y = U/V. , (27) 

si si 

and the atmosphere-related part of the Hamiltonian may be written 

H’ a = [C_ HI1 + C T (cos M HI2 + sin M HI3)] 

A mV D L 7 

1 


where 


P A 

co e 

m 


w 


(28) 


HI1 = 

UA + rA 


u rp 

HI2 = 

rA - UA 


u rp 

HI3 - 

-A V./U . 


X l 


This concludes the presentation of the basic equations for the determination of 
the optimal control and the derivation of the adjoint differential equations. 


EXTREMAL CONTROL 

In this section, applying the maximum principle, the extremal control 
will be derived for both maximization and minimization of a pay off function 
(i. e. , minimization and maximization of the Hamiltonian with respect to the 
controls) , 


Thrust Control 

Thrust Magnitude. Since the thrust is linear in the mass flow rate b, the 
expression in parentheses in equation (22) acts as a thrust switching function 

U 

SWT = — w - A . (29) 

mm v ' 

The extremal value for b is given in Table 1. 
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TABLE 1. EXTREMAL FLOW RATE 



Pay Off Maximization 

Minimization 

Value of SWT 

min 
b H 

max 
b h 

> 0 

b = 0 

b = b 

max 

= 0 

indeterminate (singular-arc) 


< 0 

b = b 

max 

b = 0 


The case of a singular arc will not be considered here. 



Figure 4. Extremal thrust 
direction and primer 
vector. 


Thrust Direction. Define the two 
vectors (Fig. 4) 

w? = (A 2 + A 2 (30) 

i ■' u rp ' 


and 

w’ = IwJ + (\/ u ) 2 ] 1/2 

A, 

then the following relations hold: 
X = w| cos <p , 

A = w5 sin <p , 

rp 1 

wJ = w’ cos T , 

A /U = w’ sin r 

X 


(31) 


(32) 


Introducing these relations into the second and third equations of equation (22) 
yields 

w = w’ (cos t cos e + sin t sin e) = w’ cos (r - e) (33) 


and 


w i = w| (cos cp cos u + sin <p sin u) = wJ cos (<p - <r) . (34) 
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The term Wj has its maximum value w| for <r - <p and its minimum -wj for 
a = n + <p. w, known as primer vector, has its maximum w’ for e = t and 
its minimum -w’ for e = ir + t. 

The extremal controls a in the vertical plane and e normal to it, 
together with the trigonometric expressions in which they appear in the 
differential equations, are given in Table 2. 


TABLE 2. EXTREMAL THRUST DIRECTION 


Case 

Pay Off 

ExtremalN^ 

Values 

Maximization 
min 
<t, e 

Minimization 
max 
O', € H 

O’ 

tt + tan -1 (A /A ) 
rp u 

tan -1 (A /A ) 
v rp u 

€ = 

7r + tan _1 |\^/ (Uwj)J 

tan-^Uw,)] 

cos e sin a = 

W(-w) 

rp 

A /w 
rp 

cos e cos cr = 

A /(-w) 
u 

A /w 
u 

sin e 

A /(-Uw) 
X 

A^/ (Uw) 

Note: The terms Wj and w are taken from equations 
(30) and (31). 



Figure 5. Deter- 
mination of ex- 
tremal bank 
angle. 


Control of the Aerodynamic Forces 

Introduce (Fig. 5) 

w 2 = [(H2) 2 + (H3) 2 ] 1/2 ; 

P = tan -1 (H3/H2) , 


(35) 


(36) 
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and the lift-dependent part of the Hamiltonian in equation (23) may be written 

H* = -C HI - C w 2 (cos (3 eos /x + sin (3 sin n). (37) 

AD L 

For positive lift coefficients C > 0, the Hamiltonian is maximal with respect 

JL 

to M for cos (j3 - /x) = -1; 1. e. , n = n + (3, minimal for n = (3. 

The lift coefficient will be confined to the range 

0 “ C t - C t • ( 38 ) 

XJ Xj 

max 

Downward accelerating aerodynamic forces may be generated by bank angles 
1(3 1 > 7r/2. 



Lift coefficients off the boundaries are deter- 
mined from 




HI - (±w 2 ) = 0. (39) 


For a drag polar of the form in equation (16) 
there follows 







(40) 


with 
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Figure 6. Determination of 
extremal lift coefficient 
from the components of 
H^ [seeeq. (37)]. 


as optimal control for C 3 = 0. 

From Figure 6 it is seen that C can 

Lt 

only be off the boundary when w 2 • HI < 0,* 
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i. e., when they have different signs. C = C for w 2 , HI 2 S 0 and 

L L 

max 

maximization of the Hamiltonian or w 2 , HI s 0 and min. H. = 0 for 

w 2 , HI ^ 0 and minimization or w 2 , HI — 0 and maximization of the 

Hamiltonian. The case C = C may also occur for w 2 • HI < 0. 

L JL 

max 

Table 3 summarizes the extremal control for the aerodynamic forces. 


TABLE 3. EXTREMAL CONTROL OF AERODYNAMIC FORCES 


Case 

Control | 

Pay Off 

Maximization 

min 

c L>f . H 

Minimization 

max 

c T , H H 

1 

•* = i 

tan” 1 (H3/H2) 

ir + tan" 1 (H3/H2) 

sin m = | 

H3/w 2 

H3/(-w 2 ) 

COS jU = 1 

H2/w 2 

H2/(-w 2 ) 

1 

C = , 






lJ w 2 , HI ^ 0 

0 

^L 

1 

| 


max 

2. | w 2 , HI & 0 

C L 

0 

1 

max 


1 

3.| w 2 , HI = 0 

singular 

singular 

4., w 2 • HI < 0 

Eqs. (40) and (41) 

Eqs. (40) and (41) 

I 

I 

with + w 2 

with - w 2 


ADJOINT EQUATIONS 

The adjoint variables X used in eqs. (20) to (28) to form the Hamiltonian 
function and in the preceding section to obtain the extremal control are given by 
the differential equations 
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'X = 


(42) 


In this section the additive terms for the inertial, thrusting, and atmospheric 
flight parts are derived. These summands of the right-hand side of the 
differential equations will be marked by the indices I, T, A. How many of 
these terms are computed and added up is controlled by logical variables from 
the main program. If all terms are used, we have 

X = \ + X T + \ * ( 43 ) 


Inertial Terms 

From equations (2l) and (42) follows: 

I “ -[ (H4 > + m rp] /r 1 


X T = U\ /r - X , 
rp, I u r 


| (H4)U - 2\ rp 1 - 2j (^ 2 (3 sin 2 A - 


A, I 


— (^ A sin X. - (HIS) cos x.) , 

U sin x \ 

b— - (A + A sm A) 

r cos 2 A v x 0 


X J 4 -— Q 6 sin A cos A , 

rp r 4 


m, I 


0 , 


'Y 1 - 0 • 
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Thrust Terms 

From eq. (22) we have for vacuum thrust 

(51) 

(52) 


U A 

e x sm € 


'\s, T b m U 


U 


m. 


U 

— b — y (±w) 
nr 


Atmospheric Terms 

Thrust Terms . From the second term in equation (22), we have 






A 

e 

3p 

*03 

3h 9h _ 

(53) 

r, 

AT = 

+ 

w 

m 

9h 

9r 5 9r “ 1 » 





A 

e 

9 P«, 

9h , 


*A. 

A, 

AT _ 

+ 

w 

m 

8 h 

TK 

(54) 





A p 
e < 

» 



•A 

m, 

AT - 

+ 

w 


"" 

• 

(55) 


There are two sets of adjoint equations for the atmospheric part, one 
for aligned inertial and atmospheric velocity vectors, which is relatively simple 
and will be treated first, and one for the general case. 

Aligned Velocity Vectors . Equations (25) to (28) yield 

HI6 = C (UA + r A ) 

D rp 

HI1 


+ C f cos m (A r - A U) - sin ju A /cos y 1 
L L u rp X J 


(56) 


HI2 
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and 


d (HI6 ) _ dC D 9Ma , v 

“au“ " C d\i + 9Ma TET lHI1 ' 


X V 

- c L ^ rp cos M - tan 2 y sin p J , 


( 57 ) 


ijHK) .c, + » (HU) 

9r D rp 9Ma dr 


+ C T f X cos p - — — tan y sin p 
u V. 


) • 


d (HI6 ) 
9r 

9 (HI6) 
9 A 


9C 


D 9Ma 


8h / 


\ 9Ma 9h 9h 

.8 (ms) 9h 
dr 9A 


(HIl) 


(58) 

(59) 

(60) 


9V. 


9V. 


W = cos Y 5 "9 F = sin 7 * 


(61) 


For the additive terms in the adjoined differential equation according to 
equation (42) 


"X = [cos y 

u, A 2m L 

- sp r • 

'X . - sin y 

rp, A 2m L 

SV. 


+ V. • 9(HI6)/9U J , 
+ V. • 9 (HI6 ) /9 r J , 


r, A 2m 
SV i 

^A , A 2m 


|f (HI6) + p • 9(HI6)/9r 


|| (HI6) + p • 9(HI6)/9A 


(62) 

(63) 

(64) 

(65) 
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( 66 ) 


*\n, A 


SpV. 



(HI6) 


General Case. Since the misalignment angles and the aerodynamic 
velocity vector are relatively complicated functions of the state variables, the 
partial derivatives of the abbreviations given in equations (23) and (24) are 
determined in the appendix. Up to equation (A-49) no adjoint variables are 
involved. The equations (A-50) to (A-56) have to be computed for each set of 
adjoint equations. 

The aerodynamic force terms in equation (23) may be written as 


• p (r, A) • (U, r, X., r, A) • H6 (U, r, X., r, A) . 

From this, the additive terms in the adjoint differential equations to U, r, and 

X. are 

i 


SpV 

* _ a 

^X, A 2m 

to r and A we have 


2 T5T H6 + v a TT 


X, A 

Finally, there is 


SV 

a 

2m 


(II (H6) + p if ) V a + 2p < H6 > 


8V 

a_ 

BX 


(67) 


( 68 ) 


m, A 



V 2 (H6) 
a 


( 69 ) 


TERMS USED IN GRADIENT METHODS 

In a general form, equations (9) to (15) are written 

X - f (x, u, t) . (70) 

The control u is eliminated in extremal field optimization methods 
by applying the calculus of variation or the maximum principle using 
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Lagrangian multipliers (see preceding sections, Extremal Control and Adjoint 
Equations) . In gradient methods, the iteration centers around an estimated 
time history for the controls. The gradient H needed in these algorithms 


u 


and the time-varying weighting matrix H -1 used in the min-H method [ 4] will 
be given in this section. 


With 


T 

X = 

(U, X., r. 

r, A 

, m, 

0 ) 



and the control vector 






u T = 

< C L> 

M, (7 

J €9 ) j 



(71) 

the matrix G 

= f 

u 

has the following form: 




§u 

§12 

gl3 

§li 

§15 




§21 

§22 

0 

§2 A 

§25 




§31 

§32 

§33 

§3 A 

S35 



G - 

0 

0 

0 

0 

0 

9 

(72a) 


0 

0 

0 

0 

0 




0 

0 

0 

G 

-1 




0 

0 

0 

0 

0 



where 








Sll = 

_Sq 

cos 

Ax a ( 

< 

cos y 

V 

Q 

0 

CO 

+ sin y & cos ix ^ 


” m 

a 8C 

L 



+ sin Ax 

a 

sin ix 

9 





§12 = - 


— C (cos u sin Ax - sin p sin y cos Ay ) 
m L a a a 
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§13 

§14 

§15 

§21 

§22 

§24 

§45 


-T 

— cos e sm or 
m 


m 


sin € cos o - , 


U 


m 


cos € cos or , 


3 C, 


Sq T / ’ ~D 

— l COS Y ■ 
Um \ 'a 3 C t 
L Li 


cos *7T- + sin y cos ji ) sin Ax 


a 




- cos Ax sin ju 
a 


_ Sq 


TT C_ (sin m sin y sin Ay + cos \x cos Ax ) 
Um L a a a 


mU 


cos e , 


U 

e 

mU 


sin e , 


§31 

§32 

§33 

§34 

§35 


Sq / , D \ 

~ l Smy a SC" - COSY a CGSM J ’ 


Sq ~ 

1 C T cos y sin \x , 
m L a 


T 

— COS € COS or 

m 


T . 

- — sin e sin o - 
m 


U 

a 

- — cos e sin a 
m 


9 


(72b) 
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From equations (22) and (23) , the following relation for is obtained: 

dC 


h t = 

u 


HI • + H2 cos p + H3 sin p 

9c l 


Sq 

m 

C T (H3 cos p - H2 sin m) 
m L 


T 

— cos € (X cos cr _ X sin a) 
m rp u 


m / X 
T / X 


(■ 


m \ U 

u 

a 

— w - X 

m m 


cos e - wj sin 


in e ^ 


which yields for H 


H = 
uu 


uu 

h ii 

^21 

0 

0 

0 


^21 

^22 

0 

0 

0 


0 

0 

h 33 

^43 

h 53 


0 ' 0 

I 

0 I 0 
I 

h 34 I h 35 

I 

h 44 i h 45 

J 

h 54 0 


(73) 


(74) 


with 


h n = - 


Sq 

m 


(HI) 


3 2 C 


D 


9C{ 


h^ = h i2 = (H3 cos p - H2 sin p) , 


h™ = — C_ (H2 cos p + H3 sin p) , 
m L 
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h 33 
h 34 = 



_T 

~m 

^43 ~ 
h 53 = 
^54 = 


e Wj 

sin e (X sin a - X cos a) , 
u rp 

U 

— — cos e (X sin a _ 
m u 

U / X 

Oi Y 

— I cos e - wi sin 
m \ U 1 


X cos a ) 
rp 


cos 

_T 

m 


9 


It is seen that thrust and aerodynamic forces are completely decoupled. Since 
the Hamiltonian is linear in the mass flow rate b, its influence may be accounted 
for by changes in the switching times t and t and the fifth row and column of 

O 1 

H can be deleted, 
uu 


For the computation of the influence functions [ 5] , the expressions 
T 

G X are needed which are similar in structure to H but are formed with 

u 

the adjoint variables for one given final condition; where as, is formed with 
the Lagrangian multipliers as a superposition of the adjoints [4] 

X = X^ + V, X T . (75) 

<t> i v l 

with v as constant multipliers for the given final conditions 4^ (X^, t^) = 0 . 


CONTROL DEPENDENT STATE SPACE CONSTRAINTS 

If the trajectory is bound to lie in a certain region of the state space, it 
is in general made up of arcs on the state space boundary and arcs off it. The 
case where the state space constraint function involves one of the controls 
directly is of special interest in atmospheric flight mechanics since both normal 
acceleration and reradiative heating constraints are of this type. 

In this section the necessary relations for applying the constraint 
optimization technique of reference [ 6] are derived for the above-mentioned 
constraints. 
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Reradiative Heating Constraint 


F 


It is assumed that a quasisteady heating equilibrium model is a good 
enough approximation to the dynamic process. Then the constraint function 
may be written in the form 

C = h (r, A) - h L j^C L , V^U, X., r, r, A)J B: 0 , (76) 

where h is a limit altitude above which the vehicle has to stay for a given 
L 

velocity and lift coefficient in order not to exceed the heating limit. On a 
boundary arc, the lift coefficient as the only control involved in the constraint 
directly has to be determined by C = 0. Then, in order to satisfy the neces- 
sary conditions of optimality, additive terms in the adjoint differential equations 
have to be evaluated (second summand) : 



From equation (76) follows 


9C 

9u 


8h i 

9C, 


0 0 0 0 


(78) 


9C _ 

9 h 

Lj 

" 9V 

a 

9V 

a 

9 V , 

i 

,8V a 

9h/9r \ 

/ 9 V 

9h/9A \ 

9X 

” 9 V 

a 

9U 

8 *i 

9r \ 

i 8r 

" 9h /9 V J 
L a / 

— 

\ 9 A 

- 9h_ /9V ) 

L a /_ 


(79) 


Together with equation (72), the second summand in equation (77) may be 
written 


8h L 

\i, o, = W K e 11 + \ + X rp 6l) 


and with the abbreviations in equation (23) 

(H7) 


_ 9C 9 V 
Sq L a 


u, Ci 


m 9V 9U 
a 



(80) 
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where 


9c 

H7 = (HI) + cos M (H2) + sin p (H3) . (81) 

8c l 

In extremal field methods, the last two terms sum up to w 2 [ see equations (35) 
to (37)] . For the other adjoints, one obtains similarly 




Sq 

~m 

9C 

9V 




•x 

x» 

Ci 

L 

9V 

a 

a 

ax, 

(H7) 

J 

(82) 




9C_ 

9V 




•x 

= 

Sq L 

a 

(H7) 

J 

(83) 

TP) 

► C 1 

- m 9 V 

a 

9r 




Sq / 

f 9C_ 

9C 

9V > 

V 


•X 

= 

L 


L a 

)• H7 , 

(84) 

r. 

C 1 

m \ 

V ah L 

" 9V 

9r 

a j 

/ 

9 *\ 


Sq i 

( 9h/9 A 

8C l 

9V \ 

3 1 - H7 

(85) 

A a 

A, 

Ci = 

In | 

V 8 V 8c l 

“ 9V 

a 

9 A J H7 * 


The derivatives 9C/9V and 9C_/9h T have to be furnished by the subroutine 
L a L L 

which generates the control on the boundary C _ (p, V ). 

J_i a 


Normal Acceleration Constraint 


The constraint equation is 

C = — p V 2 C T - n < 0 
2mg 0 a L zmax 

From 


9C 

9u 


' SpV 2 

- 0 0 0 0 

2mg 0 


( 86 ) 


(87) 
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and 


ac 

9X 


SV pC_ 
a L 

gm 0 


av 

av 

av 

/ av 

a 

a 

a 

1 

au 

8X, 

a* 

\ dr 


V q \ /8V 
a 9p \ / a 

2 p 9ry y 9A 


+ 



9p\ 

9A/ 


» 


( 88 ) 


the additive terms in equation (77) become 


•X 

u, 


C 2 


av 

- 0 * 5 S ^ V a C L W (H7) ’ 


av 

c 2 = - 0 * 5 s V g l ax. 


a 


> 


9V 

•X _ = -0.5 SV pC (H7) 

rp, C 2 a L 8r 


/9V V . \ 

= -0.5 SV pC + — — • ) 

a L\ 9r 2p 9r J 


r, C 2 


H7 


= -0.5 SV pC, 


( 


av v a . 

a + — — ^ IH7 


3 1 


"A , C 2 a L V 8 A 2p a A; 

On the boundary, die lift coefficient is determined from 


C_ = 2mg 0 n 
L u zmax 


) • 


(89) 


(90) 


Note that for C — C and C T = 0, this constraint can always be 

L Li L 

min min 

satisfied, whereas the heating constraint above may lead to a reduced set of 
entry conditions onto the constraint for feasible trajectories. This difference 
comes from the fact that the normal acceleration is generated by the control, 
while the heating constraint is a state space constraint where the control enters 
as a parameter. 
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FORTRAN SUBROUTINE 


The mathematical model described in the previous sections has been 
coded in FORTRAN for use in an extremal field algorithm. The program is 
self-explanatory once the logical parameters are described, which is the 
purpose of the following section. 

Controlling Logic 

The purpose of incorporating ten logical variables was to achieve a 
flexible and efficient subroutine. Five of the control variables are input 
parameters into the main program and remain constant throughout the iterations. 
They partly specify the case to be run. 


TABLE 4. LOGICAL CONTROL PARAMETERS 
FOR CASE SPECIFICATION 


Logical 

Variable 

State 

Result (deg) 


.True. 

Planet oblateness taken into account: altitude 

Lobl 


h = r - R 0 + Ah(A) , equation (18) 


. False. 

Spherical planet h = r - R 0 


. True. 

Planet rotation taken into account 

Lrot 

. False. 

Nonrotating planet (co^ = 0) 


. True. 

Horizontal winds W (h) taken into account 

LWI 

. False. 

No winds 



= Lrot .or. LWI 

LDVA 

.False. 

Special path for atmospheric terms to reduce 
work load in the case of aligned inertial and 
aerodynamic velocity vectors 


. True. 

Payoff maximization 

LMIMAX 

.False. 

Payoff minimization 
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The other five logical variables are stored in a two-dimensional array 
and controlled by a switching function subroutine which sets them as a function 
of the state and adjoint variables during the iteration. The first index refers to 
the trajectory section under consideration, the second to the number of the 
switching function. Their effect is given in Table 5.* 


TABLE 5. LOGICAL PARAMETERS 
SET BY SWITCHING FUNCTIONS 


Logical 

Variable 

State 

Result 


. True. 

Atmospheric terms are taken into account 

LU(J, 1) 

. False. 

Atmospheric terms are by-passed 
(e. g. , for h > 100 km) 


. True. 

Thrust is on 

LU(J, 2) 

.False. 

All thrust terms are by -passed 


• True. 

V ehicle flies on heating constraint boundary 

LU(J, 3) 

.False. 

Heating constraint not effective 


. True. 

Vehicle flies on normal acceleration boundary 

LU ( J , 4) 

.False. 

Normal acceleration constraint not effective 


. True. 

Flight with C^ 



max 

LU(J, 5) 

.False. 

C determined from X and A. 
L 


* The FOR TRAN -program is documented under separate cover. Copies are 
available from: Trajectory & Optimization Theory Branch, Aero- 
Astrodynamics Laboratory, Marshall Space Flight Center, Alabama 35812, 
Attn: Mr. J. R. Redus, S&E-Aero-GT. 
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APPENDIX 


DERIVATION OF PARTIAL DERIVATIVES NEEDED TO FORM 
THE ADJOINT DIFFERENTIAL EQUATIONS 

For nonaligned inertial and aerodynamic velocity vectors, the following 
dependencies are noted: 

V NS ( Ul W NS C h(r ' A) ]} = V NS (U > r> A) • 

V Ew{ U ’ V r> A> W EW [ h(r ’ A) ]} = V EW (Ul V r> A) 1 

V ah (V EW V NS ) ’ V ah *i’ r ' A > > 

V a (V ah’ ” V a < U > X i> r ' A) > 

X a < V ah> V a> = \ < U > *• x ,' A > • 

AX a (V EW’ V NS’ V ah' V = Ax a (U> X i' r > A) • 
p[h(r, A )] = p(r, A) , 
p ro [h(r, A)] = (r, A) , 

Ma [> , h(r, A)]= Ma (U, r, x., r, A) , 

C D (Ma, h) = C D (U, r, x r, A) . 

The equations are given in the section, Equations of Motion and Hamiltonian 
Function. From these the following partial derivatives are obtained: 

North-south aerodynamic velocity component: 


V _ = U cos x + W 
NS i NS 


(A-l) 
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9 V 


NS 


9 V, 


NS 


W = cos x i 5 “9x7 


= -U sin x. 


av aw 

NS NS 


av aw 

NS NS ah 


8r 


ah 


8A 


8h 9 A 


(A-2) 


East-west aerodynamic velocity component: 


V EW = U Sil1 X i + W EW ‘ V COS A 


(A-3) 


av av 

EW EW 

-8U - = Sln x i ' ST ' = u cos x i 


av aw av aw 

EW EW . EW EW 9h . . 

— r = — rr — - or cos A ; — — — = — — — — + o)„r sin A 

9r 8h P 8 A 9h 8A P 


(A-4) 


Radius normal aerodynamic velocity component: 


V ah = ( V EW + V NS) 2 ’ 


(A-5) 


!!£h-fv !!n + v 8V ^ /" - Av 

9U \NS 9U V EW 9U // V ah cos » 

(A-6) 


av 


ah 


9 V 


= V. 


NS 


+ V, 


9 X. \ NS 9x EW 9X 

iv i 




V . = U sin Ax , 
ah a 


(A-7 ) 


9V 


DVAHDR = 


ah 


9r 


aw. 


NS 


/aw 


NS 9h 


+ V, 


{ 


EW 


EW\ 9h 


- co cos A 


)/ 


ah , 


9V 


DVAHDL = 


ah 


av. 


9A 


= V 


NS 


av 


+ v. 


EW 


NS 9 A EW 9 A 


)/\u 


(A-8) 

(A-9) 
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Aerodynamic velocity: 


V 

a 

y^ah 

av 

a _ 

8V ah 

au 

au 

av 

a 

9V ah 

8Xi “ 

8 *i 


•2 Wo 


= cos y cos Ax 


= U cos y sin Ax 


(A-10) 


(A-ll) 


(A-12) 


3V 

___ = r/V = sin y , 
8r a a 


(A -13) 


av av . av 

a ah a 

= cos y — r — ; = cos y 

8r a 8r 8A f a 

Aerodynamic flight path angle: 

sin y = r/V , 

a a 


3 sin y . av 
DSGDU “ — — - = -~ 

V* au 

a 


au 


a sin y 


DSGDRP = 


a. = I V 
9r \ a 


3sin y . 8V 

DSGDC = — = -— 2 - 

ax. v 2 ^ ax. 

i a i 


8 sin y 


DSGDR = 


DSGDLA = 


-k = -sm y cos y 

3r 'a a 


3 sin y 

8 A ~ = - sin Y a 008 Y a 


av u 

ah 

a a 

(A -14) 


(A-15) 

y cos y 

cos Ax /V , 

a a 

a a ’ 


(A -16) 

r 

V* = COS 2 

y /V , (A -17) 

a 

a a 

' cos y 

sin Ax U/V , 

a a 

a a 


(A-18) 

av . / 

8r/ V a ’ < A -»> 

5V . / 

— /V 
a A / v a 

, (A-20) 


32 



9 


(A-21) 


cos y = V , /V 
a ah/ a 


DCGDU = 


9 cos y / 9V, 9V \ 

— = ( V — _ y ) y2 

9U \ a 9U ah 9U / a 


= sin 2 y cos Ax /V . 
a a a 


9 cos y 


DCGDRP 


DCGDC = 


q . — = -sin y cos y /V , 
9r a a a ’ 


9 cos v 

a 

9X. 


= sin 2 y sin Ax U/V , 
a a a 


(A -22) 
(A -23) 
(A -24) 


9 cos Y 9V 

dcgdr - — - -J* n* y a /v a , 


(A -25) 


9 cos y 9V , 

DCGDLA = — 5J-S. = -jf sirf y a /V a 


(A -26) 


Radius normal aerodynamic misalignment angle Ax : 

Si 


sin Ax 

a 

= < V EW 008 *i 

- sin 

NS 

X.)/V . 
i ah 

» 


9 sin Ax 

r / 9v„„ 


9V _ 

DSC DU 

a 

| f EW 

cos X. - 

NS 

9U 

V 9U 

9U 


sin x. j 


(A -27) 


ah 


9V 


- sin Ax V 


ah 


a ah 9U 


'y 2 


ah 


= -sin Ax cos Ax /V , , 

a a ah 


(A -28) 


DSC DC = 


9 sin Ax 

a 

ex. 


U(cos 2 X. + sin 2 x )V , 
i i ah 


9V 


- sin Ax V 


ah 


a ah 9x. 


l J 


ah = U COs2 Ax a /V ah’ 
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DSC DR 


DSCDLA = 


cos Ax = 


DCCDU = 


DCCDC = 


DCCDR 


8 sin A X ( 3V 8V 

__ = ( v __ cos X, - sm Xl 


9V 


- sin Ax 


ah 


a 9r 


)/’ 


ah 


9 sin Ax 

a 

9 A 


- sin Ax„ 


-( 


9v ew 

— C0S X i 


9V 

sin X. 
9A l 


9V 


ah 


)/’ 


'"a 9 A // ah ’ 

(V EW sln X i + V NS 003 X i )/V ah 
9 cos Ax 


(A-30) 


(A-31) 

(A-32) 


9 U 


- cos Ax 


( 


aV EW . , 9V NS 

9U sm x i + — COS 


»,) 


9V 


ah 


a 9 U 
9 cos Ax 


V , = sin 2 Ax o /V o , , 
ah a ah 


(A-33) 


ax. 


U(cos X. sin X. - sin x. cos X.) 


- cos Ax 


a 9x. 


ah 


-U sin Ax cos Ax /V , 
3 a. aii 


(A -34) 


8 cos A X ( 9V ew 8V ns 

57 V - 57“ sm x i + -TT- cos x i 


9V 


- COS Ax 


ah 


a 9r 


ah 


(A-35) 
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DMADR = 

9Ma _ 
8r 

/9V 

( — 

\ 8r 


(A -44) 

DMA D LA = 

8Ma ) 
8 A 

/ 8V 
^ 9 A 

- Ma l l)/ a • 

(A -45) 


Drag coefficient: 
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(A -46) 


C D (Ma, h) 


DCDDU = 


DC DDR = 


DCDDLA = 


9C 

D _ 

8C d 

9Ma 

9U 

9Ma 

9U 

8 c d_ 

8c d 

9Ma 

9r 

9Ma 

dr 

8C d 

8c d 

9Ma 

9A 

9Ma 

9A 


similar for r, X. » 


9C 


D 
9h 


^9h 
9h 9 A 


(A -47) 
(A-48) 
(A-49) 


With these expressions the partial derivatives of the composite terms of 
equation (23) may be written: 


9H0 8 003 A *a 

DHODU =— = 5,5 X 


. + A ( 
u X\ 


9 sin Ax 


a 


X\ 9U 


- sin Ax /U 
a 


)/u 

(A-50) 


Let X stand for x. > r, A ; then we have 

9H0 , 8 003 AX a , X x 8 3in AX a 

— A LVJ + 


9X 


u 9X 


U 9X 


(A-51) 


9H3 8 3in Ax a X X 8 003 Ax a 

9X u 9X “ U 9X 


(A-52) 


With X standing for U, r, x^ r, A, the partials of HI and H2 may be written 


9H1 

9X 


9H2 

9X 


9H3 

9U 


9 HO ^ „„ 8 003 y a 

_ cos y a + HO 


9H0 . . 8 Sin y s 

9X 3m y a H ° 8X 


+ X 


9 sin y 

_a 

rp 9X 
9 cos y 


- A 


rp 9X 


(A-53) 
(A -54) 


9 sin Ax 

‘ \ au ~ - \ 


( 


9 cos Ax 

a 

9U 


- cos Ax /U 
a 


)/u . 
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Finally we have as partials of H6 [see equation (24)] 


9H6 = 
9X 


aC D 

9X 


HI + C 


9H1 
n 9X 


+ C cos M 


9H2 

9X 


+ C, 


sin ix 


9H3 

9X 


(A-56) 
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